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Abstract 

In this paper, a new family of implicit compact finite difference schemes for 
computation of unsteady convection-diffusion equation with variable con- 
vection coefficient is proposed. The schemes are fourth order accurate in 
space and second or lower order accurate in time depending on the choice 
of weighted time average parameter. The proposed schemes, where trans- 
port variable and its first derivatives are carried as the unknowns, combine 
virtues of compact discretization and Pade scheme for spatial derivative. 
These schemes which are based on five point stencil with constant coeffi- 
cients, named as (5,5) Constant Coefficient 4th Order Compact [(5,5)CC- 
40C], give rise to a diagonally dominant system of equations and shows 
higher accuracy and better phase and amplitude error characteristics than 
some of the standard methods. These schemes are capable of using a grid 
aspect ratio other than unity and are unconditionally stable. They effi- 
ciently capture both transient and steady solutions of linear and nonlinear 
convection-diffusion equations with Dirichlet as well as Neumann boundary 
condition. The proposed schemes can be easily implemented and are applied 
to problems governed by incompressible Navier-Stokes equations apart from 
linear convection-diffusion equation. Results obtained are in excellent agree- 
ment with analytical and available numerical results in all cases, establishing 
efficiency and accuracy of the proposed scheme. 
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1. Introduction 



In this paper, we consider the unsteady two-dimensional (2D) convection- 
diffusion equation 

dcj) d 2 (j) d 2 (j) , . <90 ,, ,d<f) . , , , _ /n 

>rc(x,y,t)—+d(x,y,t)— = s{x,y,t), {x,y,t) e i2x(0,TJ 



or aor a|/ z <xc ay 

(1) 

for the unknown transport variable 4>(x, y, t) in a rectangular domain f2 C M 2 , 
(0,T] is the time interval, with the initial condition 

<f>(x, V, 0) = <f>o(x, y), (x, y) eCl (2) 
and boundary condition 

bi(x,y,t)<j) + b 2 (x,y,t)^ = g(x,y,t), (x,y) e dSl, t e (0,7]. (3) 

Here a is a positive constant, c(x,y,t) and d(x,y,t) are convection coeffi- 
cients, and s(x,y,t) is a forcing function which together with <p (x,y) and 
g(x,y,t) are assumed to be sufficiently smooth. b\ and 62 are arbitrary coef- 
ficients describing the boundary condition as a Dirichlet, Neumann, or Robin 
type in the boundary normal direction n. 

Equation (JT|) describes convection and diffusion of various physical prop- 
erties such as mass, momentum, heat, and energy It is encountered in many 
fields of science and engineering and often regarded as the simplified ver- 
sion of Navier-Stokes (N-S) equations. Numerical solution of the convection- 
diffusion equations plays an important role in computational fluid dynamics 
(CFD). 

The emergence and growing popularity of compact schemes have brought 
about a renewed interest towards the finite difference (FD) approach. As 
such a great deal of effort towards numerical approximation of convection- 
diffusion equation using compact FD approach can now be seen in literature. 
A compact FD scheme is one which utilizes grid points located only directly 
adjacent to the node about which the differences are taken. These schemes 
offer higher accuracy even when the grid size is small. They are able to de- 
termine the flow with information solely from the nearest neighbours. The 
major advantage of compact discretization is that it leads to a system of 
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equations with coefficient matrix having smaller band width as compared to 
non-compact schemes. For solving convection-diffusion equation, the high 
order compact (HOC) discretizations have been utilized in a number of dif- 
ferent ways and a variety of specialized techniques have been developed. The 
pioneering works on unsteady 2D convection-diffusion equation were done by 
Hirsh [l] and Ciment et al. [2] . In their works they have discussed condition- 
ally stable compact FD schemes which are fourth order accurate in space and 
second order accurate in time. In 1984, Gupta et al. 0] proposed a compact 
FD scheme for the steady convection-diffusion equation with variable coef- 
ficients. Noye and Tan |4J developed a spatially third order and temporally 
second order accurate nine point implicit HOC scheme with a large region of 
stability for 2D unsteady problem with constant coefficients. The last part of 
the previous century and early part of this century have seen significant de- 
velopment of compact algorithms for convection-diffusion equation in general 
and N-S equations in particular. In this regard the contributions of Dennis 
and Hudson @], Spotz and Carey @, [J, Li et al. Q, Zhang Q, Kalita 



et al. [10J, Dehghan [TjJ deserves special mention. In the last decade or 
so HOC schemes for convection-diffusion equation have also been developed 
for nonuniform grids 



12, 13, 14 . Of late a lot of efforts have been made to 



increase computational efficiency of HOC schemes. This has led to the devel- 
opment and implementation of HOC-ADI algorithms 15|, Ll6|, Ll7|, ll8|] - These 
HOC-ADI methods combine computational efficiency of ADI approach and 
high-order accuracy of HOC schemes. Besides, various multigrid and conver- 
gence acceleration techniques have been investigated and proposed to solve 
convection-diffusion equations discretized by the fourth order compact finite 
difference schemes 



191, 120|, I2JJ 



The aim of this paper is to introduce a new family of implicit compact FD 
formulation for solving unsteady convection-diffusion equation with variable 
convection coefficients. The developed schemes are fourth order accurate in 
space. One of the schemes also possess second order accuracy in time and rest 
are first order accurate. The approach involves discretizing the convection- 
diffusion equation using not only the nodal values of the unknown transport 
variable but also the values of its first derivatives at a few selected nodal 
points. Inturn first derivatives are discretized by using Pade approximation. 
The new family of schemes are established on a rectangular domain using only 
the four nearest neighboring values of <fi. Although developed for convection- 
diffusion equation with variable coefficients the schemes are accessible in the 
form of a diagonally dominant system with sparse constant coefficient matrix 
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and are named as (5,5) Constant Coefficient 4th Order Compact [(5,5)CC- 
40C] schemes. This significantly decreases computational complexity. In 
the process we also remove the usual restriction on HOC schemes of having 
to use a grid aspect ratio of unity. At first we develop a new HOC scheme 
for the ID steady convection-diffusion equation to carry out wave number 
analysis and compare with other HOC schemes. Subsequently the family 
of implicit 40C schemes are proposed and are shown to be unconditionally 
stable in von Neumann sense. To validate the accuracy and efficiency of the 
newly derived compact schemes numerical experiments are performed on five 
test problems. They are: (i) Taylor's vortex problem, (ii) The convection- 
diffusion of Gaussian pulse, (iii) Analytic solution of N-S equations with a 
source term, (iv) Doubly-periodic double shear layer problem, and (v) Lid 
driven cavity flow. Excellent comparison can be seen in all the cases. 

The rest of this paper is organized into three sections. In section 2 we 
outline the (5,5)CC-40C family of schemes and its properties. In section 3 
we present numerical results and finally in section 4 concluding remarks are 
offered. 

2. The (5,5)CC-40C schemes 

2.1. Introduction of the new approach for steady ID convection- diffusion 
equation 

In order to describe basic ideas of our approach we start with elementary 
ID convection-diffusion equation 

- <Pxx + c(x)4> x = s(x). (4) 

A fourth order compact approximation of (j3J) is obtained by considering the 
discretization 

^ = 2<*&-<$ <e Bi + O(/i*) (5) 

where 8 X and b 2 x are the second order central difference operators for first and 
second order derivatives respectively, h is the mesh length along x- direction 
and <f>i denote the approximate value of 4>{ x i) at a typical grid point X{. Using 
(j3J) in (j3J) and retaining the leading truncation errors we obtain a fourth order 
compact discretization 

- + ician) + 5 x )<j> Xi = Sl + 0(h 4 ). (6) 
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The main advantage of the discretized equation is that it leads to a 
diagonally dominant tridiagonal system with constant coefficient in transport 
variable cf) and as such quite efficient for numerical computation. But the 
main disadvantage is that the derivative (f) Xi is required to be approximated 
upto fourth order of accuracy. This is separately achieved by using standard 
Pade discretization 

4>Ti = M-^M + o{h*). (7) 

o 

The classical truncation error analysis shows that the relation ([6]) along with 
(J7D constitute a compact fourth order discretization of the differential equa- 
tion (J3J). The idea of using derivative of flow variable in addition to the 
variable itself in the process of discretization, although new for convection- 
diffusion equation, has its roots in the work of Gupta and Manohar 22J on 



biharmonic equation. Later it was formulated by Stephenson [23 



2.2. Modified wave number analysis 

To examine in detail the characteristics of the fourth order compact 
scheme developed above, a modified wave number analysis of equation (TJJ 
is carried out. Note that truncation error does not represent all the charac- 
teristics of a finite difference scheme and the modified wave number analysis 
provides a tool to compare the resolution of finite difference operator with its 
analytical counterpart. It allows one to assess how well different frequency 
components of a harmonic function in a periodic domain are represented by 



a discretization scheme [161 . Il8j . Wave number analysis also provide an op- 
portunity to compare this new philosophy vis-a-vis with some of the well 
established finite difference algorithms. 

Let us consider the trial function e lKX (I = y/—l) on a periodic domain 
where k is the wave number. The exact characteristic of the differential 
equation fll]) is 

^Exact = K 2 + ICK. (8) 

We determine the characteristic of scheme ([6]) by replacing the difference 
operators with the corresponding wave number and using the equation (J7]). 
The characteristic of the newly proposed fourth order compact scheme is 
given by 

5 — 4 cos(k/i) — cos 2 (/t/i) 3sin(/c/i) 
^cc-aoc = h 2 {2 + cos{Kh)) + C /i(2 + cos(^))' ^ ' 
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For the sake of comparison the characteristics of some of the well known 
schemes are presented below. 

The characteristic of the second order central difference (CD) scheme 

(-61 + c5 x )fr = Si (10) 

for convection-diffusion equation (J3J) is 

Acd = Ai + Ic\ 2 (11) 

2 — 2cos(k/i) sin(/t/i) 
where A x = — , A 2 



The characteristic of fourth order HOC based schemes 



(-ai<^ + c^ x )0j = (1 + a 2 6l + a 3 6 x )si (12) 

is 

a 1 X 1 + Ic\ 2 , , 

W " (l-a^ + ZasV (13) 

Pe 2 1 — «i /i 2 1 — ai 
where a x = 1 + — — , a 2 = 5 h — , a 3 = 

12 (TO c 

and Pe = c/i. Pe is known as the cell Reynolds number. 

The characteristic of the fourth order Pade (PDE) approximation [l 

- 4>xx t + C(j) Xi = Si (14) 

12 

<j>xx. i+1 + lO0 raj + <$> XXi _ x = 75(^-2^+^) 



^2 ^+1 

3 

s 1 



12(1 — cos(k/i)) 3sin(/c/i) 
P1)i? = /i 2 (2 + cos(^)) + C /i(2 + cos(^))' (15) 



IS 



Finally the characteristic of the RHOC scheme pi 

(-/3 X 5 2 + cSJfa = (1 + /3 2 5 2 + 3 6 x ) Si (16) 

is 

/3iAi + /cA 2 , , 

W " (i-AAO + m' (17) 
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f±# + f c^O f i^i c ^o 



where ft = p 12 . ^ , & = < , # 
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1 — -I- — 1.2 

1 6+36 lg c = 10 C = 0. 

Figure [U shows the non-dimensionalised real and imaginary parts of char- 
acteristics (As) as function of nh. We non-dimensionalise real and imaginary 
parts of As by multiplying with h 2 and h/c respectively. Note that character- 
istics of the schemes ffT2l) and (JIB]) depend on Pe. Here we have considered 
two different cell Reynolds numbers Pe = 0.1 and Pe = 100. For the case 
Pe = 0.1, the real part of the characteristic shown in figure (Ha), the HOC, 
PDE and RHOC schemes are indistinguishable concluding identical dissi- 
pation errors whereas the second order accurate CD scheme shows much 
larger dissipation error. The figure clearly depicts that in comparison to 
all these schemes resolution of characteristic by the present compact scheme 
is much better and possess smaller dissipative error. The non-dimensional 
imaginary part of the characteristic corresponding to Pe = 0.1 presented 
in figure QJb) shows that dispersive error of the present scheme is same as 
that of PDE scheme. This is also evident from equations (Q and ( I15p . Both 
theses scheme show much smaller dispersive error when compared to HOC, 
RHOC and CD schemes. Note that at Pe = 0.1 the imaginary part of the 
characteristic for HOC and RHOC schemes are basically same. As the cell 
Reynolds number is increased to Pe = 100, the dissipation error for the HOC 
and RHOC scheme increases dramatically as shown in figure He), while the 
newly developed present scheme as also the PDE and the CD schemes do 
not alter their nondissipative properties. Again one can see the superiority 
of the present scheme. In figure [TJ(d) we compare the dispersive errors of the 
schemes at Pe = 100. One can clearly see a significant overshoot produced 
by HOC scheme. Here the present scheme produces the same resolution as 
that of the RHOC and the PDE schemes. The above discussion suggests 
that the newly developed compact scheme, which do not alter its character- 
istics with cell Reynolds number, possess resolution properties better than 
the schemes discussed here. Thus the present algorithm may be best suited 
for high Reynolds number flow computation. 

2.3. (5,5)CC~40C scheme for the 2D unsteady case 

To begin with, we briefly discuss the development of HOC formulation 
for the steady state form of equation (TT]), which is obtained when c, d, s and 
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are independent of t. Under these conditions, equation (jTJ becomes 
d 2 <f) d 2 <f) x 9 ^ , j/ ^ i \ 

- w - w + c{x ' y) ^ + d{x ^ = s{x > y) - (18) 

In order to obtain a compact spatially fourth order accurate discretization 
we consider the following approximations for second order space derivatives 
appearing in equation ( 1181) 

<t> XXi . = 28l<t> i .-8 x( t> Xi .+0{h A ), (19) 

^o„-^„ • W). (20) 

Here h and k are the mesh lengths along x- and y- directions respectively and 
(fiij denote the approximate value of <p(xi, yj) at a typical grid point (xj, yj). 
The above approximations yield an 0(h A , fc 4 ) approximation for equation 
f|T8|) on a five point stencil as 

- 2^0 SJ - 25^. + (<\,. + Cij)(f) Xi j + (5 y + dij)^. = sij. (21) 

Compatible fourth order approximations for space derivatives <p Xi and 4> yi 
appearing in equation (I2T!) may be obtained using 

h 2 

<f>xij = ($x<f>i,j ~ y^^J, (22) 

and 

k 2 

respectively. Note that compared to standard HOC formulation |l0|, we are 
not required to assume any smoothness condition on the convection coef- 
ficients c, d and forcing function s. The equation ( 1211) can be viewed as 
a banded system with only five non zero diagonals; of course drawback of 
requiring to approximate (f) x . . and 4> Vi . separately using fT22|) and ff23l) re- 
spectively remain. 

For unsteady case, the equation with variable coefficients will be similar 
to equation ( fl8|) . but the coefficients c and d are functions of x, y and t, and 

the expression on the right hand side becomes s(x,y,t) — a ~Q^- We intend 

to discretize time derivative as accurately as possible and obtain a stable 
numerical scheme. Introducing weighted time average parameter i such that 
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t L — (1 — i)ti + iti for ^ i ^ 1, where (n) denote the n-th time 
level, we obtain a family of integrators; for example, forward Euler for l — 0, 
backward Euler for i = 1 and Crank-Nicholson for t = 0.5. Consequently the 
finite difference schemes for equation ([T]) is 



[a - 2t5t(8l + 5 2 y )]<j>^ = [a + 2(1 - i)5t{% + %))<f>W 

+ (1 - L)St[-(S, + cg^W _ (<Sy + + s g)] 

+ ^[-(.y. + c g +i) )0£: i} - (5, + + s 5 +1) i- 

(24) 



The accuracy of the schemes are 0(5t s , h A , k 4 ), with s ^ 2. The second order 
accuracy in time is obtained for i = 0.5. All the schemes arising in this way 
are implicit and in the following subsection it will be shown that this family 
of schemes is unconditionally stable for 0.5 ^ i ^ 1. For all values of i in 
this range, except i = 1, the difference stencil requires five points in both 
n-th and (n + l)-th time levels resulting in what we call (5, 5) scheme. The 
compact stencil emerging in this way has been illustrated in figure [2j 

2.4- Stability analysis 

In this section we carry out von Neumann linear stability analysis of the 
finite difference scheme (I24p by assuming the convective coefficients c and d 
to be constants and forcing function s in equation ([1]) to be zero. 

Theorem 1 (Stability). The finite difference scheme (Efy is uncondition- 
ally stable, in the von Neumann sense, for 0.5 ^ i ^ 1. 

Proof. Let <j$ = b^e Idxl e Idvj where is the amplitude at time level n, 
and 9 X = 2irh/Ai, 9 y = 2-nkj A 2 are the phase angles with wavelengths Ai 
and A2 respectively, then from relations (T221 and fT23l) we get 




3 sin 6 X 



,(«) 



(25) 



h(2 + cos 6 X ) 



and 
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6 



(26) 
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Also 



(Si + 5 2 y )^ = (|(cos0 s - 1) + liccedy - (27) 



(28) 



Using relations ( 125|) - (!29|) in the scheme fl2^|) we obtain 



. cos 6^ - 1 cos0 y - 1 
a - Atdt ( — h 



h? k 2 

o + 4(l-t)«( ^ 2 
sin ^ 



i.3 



cos 6^. — 1 cos Oy — I 



k 2 



+(1 - t)(Jt 



sin 6* r 



- Ic 



3 sin 0,, / sin 6L 

+ I — ^ - Id 



h(2 + cos8 x ) V 



(29) 



3 sin 9,, 



3 sin Q~ 



h J h(2 + cos 9 X ) 

If G is the amplification factor then 



+ 



sin 6L 



Id 



k(2 + cos 9 y )_ 
3 sin 0„ 



fc(2 + cos ^ 



(n+l) 

i,3 



1 + 



(l-t)<5f 
a 



G 



cos 2 63; +4 cos 9 x -5 I cos 2 6^+4 cos 8 y 
h 2 (2+cos6 x ) ~ t ~ fc 2 (2+cos6» H ) 



-5 



_ t( _3sin0 2 _ i j 3 sin 6> y 

1 °/i(2+cose a; ) U A:(2+cos(9 H ) 



cos 2 6^+4 cos 9 x -5 , cos 2 6^+4 cos 6> y -5 
h 2 (2+cos 9 X ) 1 fc 2 (2+cos H ) 



J C 



3 sin 6 



3 sin B v 



h(2+cos6 x ) 1 fc(2+cos9y) 



1 + {1-l)(A + IB) 
1 - + 



where 



and 



.4 



St ( cos 2 9 X + 4 cos 6^. — 5 cos 2 0„ + 4 cos 0„ — 5 



a \ h 2 (2 + cos 6* 2 



+ 



<ft / 3sin^ 
B — — c— — + d- 



k 2 {2 + cos 9 y ) 
3 sin 9„ 



a \ h(2 + cos 9 X ) k(2 + cos 9 y 
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Therefore 

(i + A-^) 2 + (i-Q 2 i? 2 

11 (l-^) 2 + i 2 5 2 

For stability we require |G| ^ 1 which gives 

2A + (A 2 + B 2 ){\ - 2l) ^ 0. 

cos^ ~\~ 4 cos $ — 5 

Since ^ V 9 G [0, 2tt1 we note that A < 0. Thus the 

2 + cosfl L J 

above inequality is clearly satisfied if 

(1 - 2i) < =>- i ^ 0.5. 

This completes the proof. 

2.5. Solution of algebraic system 

Let us now discuss the solution of algebraic systems associated with the 
finite difference approximation ( )24|) . Introducing the notations 

$ = (01,1, 01,2, 0m,n) T , §x = (0x1,1, 0x1,2, ■••> 0x m ,n) T and $j, = (0^, yi 2 , (j)y m J T 

the resulting system of equations in matrix form can be written as 

Mi$ (n+1) = $W ) $W ) (30) 

For a grid of size mxn, the matrix Mi has the dimension mn. For the scheme 
developed here Mi is a diagonally dominant constant banded matrix with five 
non zero diagonals. Also $i n) , $ (n+1) , $i n+1) , $j, n+1) are all mn 
component vectors. At any time step, once $^ n ^ has been approximated $i n ' ) , 
can be obtained by solving tridiagonal systems 

M 2 $i n) = F 2 ($ (n) ) (31) 

M 3 $< n) = F 3 ($ (n) ) (32) 

respectively. The equations f l3~T]) and ( )32|) are the corresponding matrix forms 
of the relations f )22|) and (!23|) . System ( 13T]) and ( 132|) being tridiagonal systems 
they can be solved using highly efficient numerical algorithms. Thus the main 
objective is to solve the system ( 1301) . thereby evaluating unknown vector 
Note that difficulty arises due the presence of (n + l)-th time level 
gradients of $ on the right hand side of equation ( 130]) as those quantities 
will be available only after solving for transport variable at the (n + l)-th 
time level. To overcome this difficulty we adopt correcting to convergence 
approach. The entire strategy can be summarised in the following algorithm: 
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1. Begin with $ (n) . 

2. Obtain $i n) and $J n) using ([SID and (122) respectively. 

3. Take ^ = *W ^ = *W ^ = *< n) . 

4. Correct to $( n+1 ) using (HQ}. 

new; " ' 

5. Correct to using ( 13T|) and (|32|) respectively. 

6. If ||<l>( n+1 ) - $( n+1 )|| < e then $( n +!) = $( n+1 ). 

1 1 new old ' ' new 

7- = ^ = = ^ goto step 4. 

Direct solution of any of the above linear system is impractical because of 
huge size of the coefficient matrix. We solve the system using the biconju- 
gate gradient stabilized (BiCGStab) j24| method without preconditioning, 
where, due to the compact nature of grid, it is easy to implement matrix 
vector multiplication Mi$ without the need of storing any of the entries of 
the matrix Mi. The convergence criterion for BiCGStab iteration based on 
norm of residual depend on grid size and problem under consideration but 
a typical choice may be 10 -10 . Similarly, a typical stopping criterion for the 
inner iteration can be set at e = 10 -12 . We carry out our computations 
on a Pentium Dual-Core processor based PC with 4 GB RAM using double 
precision floating point arithmetic. 



3. Numerical Examples 

We conduct numerical tests to examine and compare the accuracy and 
efficiency of the present fourth order accurate scheme. We use the pro- 
posed scheme and obtain numerical results for five different problems of 
varying complexity. The first two of these five problems are concerned with 
convection-diffusion equation and the other three solves Navier- Stokes equa- 
tions in two dimensions. For all the problems discussed here we consider 
uniform mesh (h = k) and take i = 0.5 to ascertain second order temporal 
accuracy. 

3.1. Problem 1: Taylors vortex problem 

This is a pure convection problem in the unit square domain [0, 1] x [0, 1] 
with the coefficients a = 1 and c = d = 0. The exact solution of this test 
problem [15|, ll7|, ll8| is 



(j)(x,y,t) = e n sin(7nr) sm(iry). (33) 
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In order to compute using discretized systems (1241) along with (1221) and 
the required initial and Dirichlet boundary conditions for 0, <p x and <p y can 
be directly derived from (133]) . We consider three different mesh sizes 11 x 11, 
21 x 21 and 41 x 41 and compute errors with respect to the exact solution 
using L\, L2 and norms. The computed results at time t = 0.25 and 
t = 0.5 have been presented in tabled) A spatial order of accuracy close to 
four, for all cases, can be seen in this table. We estimate order of accuracy 
using the formula lri2{Errl / Err2) where Errl and Err2 are the errors with 
the grid sizes h and h/2 respectively. Temporal second order accuracy of the 
scheme is verified in the table [2] where we compute using three different time 
steps St = 0.01, St = 0.005 and St = 0.0025 keeping h unchanged. The grid 
convergence along the line x = 0.5 have been presented in figure |3fa) which 
establishes grid independence of the scheme. In figure E£b) we present time 
evolution of Li- norm error for three different grid sizes from which stability 
and effectiveness of the scheme can be gauged. 

3.2. Problem 2: Convection- diffusion of Gaussian pulse 

To further study the effectiveness and validity of the new scheme we apply 
it to an unste ady problem concerning the convection-diffusion of a Gaussian 



pulse 15|, |l6|, ll8j in the square region [0,2] x [0, 2] with analytical solution 

(ax — ct — 0.5a) 2 (ay — dt — 0.5a) 



(f)(x,y,t) = Af + 1 ex P 



a(4t + 1) a(4t + 1) 



(34) 



The initial condition <p(x,y,0), taken from (jM|) . is a Gaussian pulse cen- 
tered at (0.5, 0.5) with pulse height 1. Dirichlet boundary conditions are 
also directly obtained from the analytical solution. In this study we con- 
sider two different combinations of convection coefficients c = d = 80 and 
c = d = 10000 with fixed value of a = 100. For the first case we compute 
using three different grids and present errors calculated using three different 
norms at time t = 0.5 and t = 1.0 in the table |3j This table again verifies the 
analytical fourth order spatial convergence of the scheme. Using identical 
combination of the coefficients we estimate the temporal order of accuracy 
of the scheme in table |H In this table we have used three different time 
steps St = 0.1, St = 0.05, St = 0.025 keeping spatial grid spacing fixed at 
h = k = 0.05. An order of accuracy around two can be seen in all the cases. 
The temporal decaying nature of error computed using three different norms 
for cell Reynolds number Pe = 2 has been shown in figure HI As expected 
all the norms show similar type of decaying character with time. 
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The second combination of coefficients viz. a = 100, c = d = 10000 has 

In figure EJ^a) we compare the 



been motivated by the works of [15j, [ly, 118 



analytical solution and the computed solution in the region 1.2 ^ x, y ^ 1.8. 
This figure clearly shows that our computed solution is almost indistinguish- 
able from the exact one. For this computation we have kept cell Reynolds 
number fixed at Pe = 200 and have used 5t = 2.5 x 10~ 5 . Figures E^b)-(d) 



obtained from [18( present similar comparisons for the schemes given by You 



16j . Karaa and Zhang 15 1 and Tian 181 respectively. It is seen that although 



schemes given by You 16| and Tian 18[ captures the analytical solution ac- 
curately but the scheme of Karaa and Zhang 15j fails to obtain the desired 
accuracy mainly due to the enhanced numerical dissipation. Note that in 
section 12.21 we have established enhanced phase error and hence numerical 



dissipation associated with the HOC type schemes discussed in [15 



3.3. Problem 3: Analytic solution of N-S equations with a source term 

Next we consider the Navier-Stokes system of equations with a source 
term given by 



duo 1 , 

— + UU X + VU y - —\LO xx 



dt 



Re 



+ CO 



mi i 



^{x 2 + y 2 + A)e- 



t 

Re 



(35) 



- (4>xx + ^yy) = CO. (36) 

This problem considered in 25j, |26J, is defined on the square domain [0, 1] x 
[0, 1] and has analytical solution ip = (x 2 + y 2 ) 2 e~R^ , to = — 16(x 2 + y 2 )e~^, 
which decays rapidly with time. The presence of u = ip y and v = —ip x 
introduces non-linearity into the system. This problem is studied to discusses 
applicability and effectiveness of the newly developed formulation in tackling 
coupled non-linear N-S system. We begin by observing that formulation ( |24|) 
can be directly applied to equation ( 135]) where to plays the role of transport 
variable. On the other hand for equation (|36j) we use steady formulation 
(12~TT) with Cij = dij — V The necessity to compute ip x and ip y using 
Pade approximations ( |22l) and (1231 alleviates the need to compute u and v 
separately although one still has to do some additional computations for co x 
and co y . At each time step we solve all systems concurrently in tandem, using 
the strategy outlined in section I2.5[ where we impose convergence criterion 
on ip. Initial and boundary conditions necessary for computation are derived 
from analytical expressions. Table [5] and M display ip and to errors respectively 
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as measured with respect to three different norms. The computations have 
been carried out by taking Re = 1. We present results for three different 
time levels and see that the convergence rate is around four for both the ip 
and uj. From the data presented in table |5] and [6] it is noted that our Ly , 



norm errors for ip and u are much better than the errors reported in [26 



At time t = 1.0, our ip and u errors are respectively about one-tenth and 



one-hundredth of the errors reported in [26 



3.4- Problem 4'- Doubly-periodic double shear layer problem 

The doubly periodic double shear layer flow problem has been used by 



authors [271 . |28| to test the behavior of discretization for non-stationary sit- 
uations where steep gradients are involved. The shear layers are perturbed 
slightly at the initial time, which causes it to roll up in time to large vortical 
structures. The unsteady incompressible N-S equations 

du 1 

— +UU X +VUy - — (UJ XX +COyy) = (37) 

ot tie 



- (4>xx + ^yy) = U. (38) 

are solved to obtain flow field. The initial conditions for u and v are given 
by 

^,0) = (tanhto-x/2)/,] K0 <y< * 

V 1 |tanh[(3vr/2 - y)/p] if tt ^ y ^ 2tt 
v(x,y,0) = crsin(x) 

from which expressions for ip and u are obtained. We consider shear layer 
width parameter p = tt/15 and size of perturbation o = 0.05 along with 
Re = 10000 which is an excellent case to verify the precision and stability of 
the numerical scheme at high wavenumbers. Periodic boundary conditions 
are imposed along both x and y directions. 

The calculated values of the horizontal velocity u along the vertical cen- 
terline and the vertical velocity v along the horizontal centerline at time 
t = 6.0 and t = 10.0 are presented in figures [6(a) and E(b) for the grids 
65 x 65, 129 x 129, and 257 x 257. These figures verify the grid independence 
of the numerical results obtained using this newly developed algorithm. We 
see that all the sharp velocity gradients have been captured accurately by the 
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grid of size 129 x 129 and is enough for accurate resolution of the flow. For 
this problem, in the absence of analytical solutions, we estimate the perceived 



order of convergence p by using the formula [32 



1*3-^1 II K-h\ 



||^3-*2|| K- hi 

Here h\, hi, h?, are mesh sizes in three different uniform grids used to compute 
solutions vectors \l/2, ^3 respectively. In table [7] we present differences 
in the horizontal velocity fields computed using three different grids 65 x 65, 
129 x 129, and 257 x 257. We calculate differences using L\ and L 2 norms 
at time t = 1.0, t = 5.0, and t = 10.0 and estimate perceived order of 
convergence by using f|40|) . An identical estimation has also been carried 
out for the vertical velocity field in table |HJ From these tables we see that 
at time t — 1.0 the perceived order of convergence compares well with the 
analytical order of convergence, but with time progressing there is a drop in 
the perceived order of convergence for both u and v which may be attributed 
to the increasing complexity of the flow field due to severe roll up of the 
shear layer. Finally stream function and vorticity fields at time t = 6.0 and 
t = 10.0 computed using 257 x 257 grid is presented in figure [3 It is seen 
that a very good resolution is obtained by using our scheme. The present 
results are smooth and the integration is stable which is also sustained by the 
figure [6j These results show that the present method is efficient for solving 
problems with steep gradients. 

3.5. Problem 5: Lid-driven cavity flow 

We turn now to the problem of 2D lid-driven cavity, which has been used 
as benchmark test problem by many authors. This problem is of great scien- 
tific interest because it displays almost all of the fluid mechanical phenomena 
for incompressible viscous flows in the simplest of geometric settings. Over 
the years, this problem has become one of the most frequently used bench- 
mark problem for the assessment of numerical methods, particularly for the 
steady state solution of the incompressible fluid flows governed by the N-S 
equations (I3"7]) - (13"8]) . In particular, we compare our results to the steady state 



results of Ghia et al. [29(, Bruneau and Jouron [30J, and Gupta and Kalita 
3jJ. The cavity is defined as the square ^ x, y ^ 1. The fluid motion is 
generated by the sliding motion of the top wall of the cavity in its own plane 
from left to right. Boundary conditions on the top wall (y = 1) are given as 
u — 1, v — 0. All other walls of the cavity are at rest. 
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We obtain the steady state solutions using a time-marching strategy. The 
steady state was assumed to reach when the maximum -0-error between two 
successive outer temporal iteration steps was smaller than 1.0 x 10~ 10 . Nu- 
merical solutions for the driven-cavity flow are obtained taking Reynolds 
number Re = 1000, 3200, 5000, and 7500. The Reynolds number for the 
lid-driven cavity problem is defined as Re = UL/u, where U — 1 is the char- 
acteristic velocity, L — 1 is the characteristic length and v is the kinematic 
viscosity. In the present computation we consider an uniform grid spacing 

with mesh size h = k = — for Re = 1000 and h = k = for the other 

64 128 
three cases. The time increment for all the cases is equal to 5t = 0.01. 

In figure El we present comparisons of our steady state results of the 

horizontal velocities on the vertical centerline and the vertical velocities on 

the horizontal centerline of the square cavity for Reynolds numbers 1000 to 



7500 with the data from Ghia et al. [29J. In each case, velocity profiles 
obtained by the proposed method on relatively coarser grids match very well 
with results given by Ghia et al. Note that for Re = 1000, 5000 and 7500 we 
have used mesh spacing twice that of Ghia et al. 

Figure [9] depicts streamline contours for Reynolds numbers 1000, 3200, 
5000, and 7500. In these graphs, the typical separations and secondary vor- 
tices at the bottom corners of the cavity as well as at the top left can be seen. 
For Re = 7500 a tertiary vortex can also be seen at the bottom right corner. 
The vorticity contours corresponding to the above mentioned Re values can 
be found in figure [101 In this figure we have clearly depicted vorticity con- 
tours in the vicinity of the solid walls as well as at the centre of the cavity. 
The stream lines and vorticity contours obtained here are in good agreement 



with the established results available in literature [28] , I29L |30j, |31j, thereby 
confirming that our formulation yields qualitatively accurate solutions. 

To further validate the present method we provide quantitative data com- 
parison of our solutions in table [91 Here we present our steady state data 
for the primary, secondary and tertiary vortex centres and compare them 



with the well established work of Ghia et al. [29], Bruneau and Jouron [30 



and Gupta and Kalita [31|. A very close comparison can be seen except 
possibly for the strength of the tertiary vortices which may be attributed to 
the coarser nature of the grid being used. It is clear from the table that the 
results of the present numerical method are reliable and the algorithm can 
be used to solve unsteady viscous incompressible flows. 
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4. Conclusion 

A new family of implicit spatially fourth order accurate compact finite 
difference schemes have been proposed with weighted time discretization for 
unsteady convection-diffusion equation with variable convection coefficients. 
The schemes are second or lower order accurate in time according as time 
discretization parameter i = 0.5 or otherwise. Linear stability analysis shows 
that for 0.5 ^ i < 1 the schemes are unconditionally stable. The novelty of 
the family of schemes lies in the fact that the second order spatial derivatives 
of the transport variable are discretized using not only the nodal values of the 
transport variable but also its first derivatives leading to fourth order approx- 
imation. Pade approximation used only for first derivatives ensures overall 
fourth order accuracy and diagonal dominance of the resultant constant coef- 
ficients algebraic systems. Modified wave number analysis confirms that the 
family of schemes enjoys better resolution properties and hence lesser dissipa- 
tion error when compared to other compact schemes. The newly developed 
family of schemes have been used to solve Navier-Stokes system also and 
show excellent performance in terms of computational accuracy, efficiency, 
and stability with smaller number of grid points. Five different numerical 
experiments performed to demonstrate high accuracy and efficiency of the 
present method produce excellent comparison. The extension of the proposed 
discretization procedure to irregular domain and three dimensional cases are 
straightforward. 
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Table 1: Problem 1: L\-, L2-, L^- norm error and spatial order of convergence with 
St = h 2 = k 2 

Time 11 x 11 Order 21 x 21 Order 41 x 41 

t = 0.25 U 3.635xl0~ 5 3~81 2.598xl0~ 6 3^91 1. 730 xl(T 7 

L 2 4.690xl0~ 5 3.84 3.274xl0~ 6 3.92 2.156xl0- 7 

Loo 9.758xl0~ 5 3.87 6.676xl0- 6 3.94 4.354xl0- 7 

t = 0.5 Li 4.277xl0~ 7 3.66 3.384xl0~ 8 3.84 2.367xl0~ 9 

L 2 5.520xl0~ 7 3.69 4.265xl0~ 8 3.85 2.951xl0~ 9 



1.150xl0- 6 3.72 8.702xl0- 8 3.87 5.960xlQ- 



9 
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Table 2: Problem 1: L\-, L2-, L^- norm error and temporal order of convergence with 
h = k = 0.05 

Time St = 0.01 Order St = 0.005 Order St = 0.0025 

t = 0.25 U 4.119xl0~ 5 L98 1.033xl0" 5 L99 2.598xl0~ 6 

L 2 5.189xl0~ 5 1.99 1.302xl0- 5 1.99 3.274xl0- 6 

Loo 1.058xl0~ 4 2.00 2.654xl0~ 5 1.99 6.676xl0~ 6 

t = 0.5 Li 5.323xl0~ 7 1.99 1.343xl0" 7 1.99 3.384xl0~ 8 

L 2 6.708xl0~ 7 1.99 1.693xl0~ 7 1.99 4.265xl0~ 8 

1.369xl0- 6 1.99 3.454xlQ- 7 1.99 8.702xlQ- 8 
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Table 3: Problem 2: L\-, L 2 -, L x - norm error and spatial order of convergence. Here 
a = 100, c = d = 80 with 6t = h 2 = k 2 . 

Time 21 x 21 Order 41 x 41 Order 81 x 81 

t = 0.5 T x 9.902xl0~ 4 4~76 3.656xl0~ 5 4T0 2.135xl0- 6 

L 2 3.052xl0" 3 4.24 1.617xl0~ 4 4.11 9.394xl0~ 6 

Loo 2.836xl(T 2 3.89 1.915xl0~ 3 4.11 1.107xl0~ 4 

t = 1.0 Li 5.363xl0~ 4 4.75 1.992xl0~ 5 4.05 1.201xl0~ 6 
L 2 1.297xl0- 3 4.23 6.904xl0~ 5 4.06 4.134xl0~ 6 
1.036xl0~ 2 4.06 6.200xl0~ 4 4.00 3.865xl0~ 5 
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Table 4: Problem 2: L x -, L 2 -, Loo~ norm error and temporal order of convergence. Here 
a = 100, c = d = 80 with h = k = 0.05. 

Time 5t = 0.1 Order g = 0.05 Order 5t = 0.025 

t = 0.5 Li 3.608xl0" 3 L84 1.006xl0" 3 T96 2.592xl0~ 4 

L 2 1.507xl0 -2 1.79 4.361xl0- 3 1.93 1.147xl0 -3 

L^ 1.551xl0 _1 2.00 4.599xl0~ 2 1.90 1.233xl0 -2 

t = 1.0 Li 3.736xl0~ 3 2.00 9.291xl0~ 4 1.98 2.354xl0~ 4 
L 2 1.126xl0 -2 1.86 3.102xl0 -3 1.97 7.928xl0~ 4 
8.383xl0~ 2 1.75 2.494xl0~ 2 2.02 6.148xl0~ 3 
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Tabic 5: Problem 3: L\-, L2-, L^- norm error in ip and spatial order of convergence with 
5t = h 2 = k 2 



Time 



11 x 11 



Order 



21 x 21 



Order 



41 x 41 



£ = 0.5 Li 6.403xl0" 8 
L 2 8.213xl0- 8 
Loo 1.661xl0~ 7 



4.610xl0" 9 
5.769xl(T 9 
1.154xl0~ 8 

2.838xl0~ 9 
3.551 xl0~ 9 
7.102xl0- 9 



3.082xl0" 10 
3.814xl0- 10 
7.565xl0- 10 

1.899xl0~ 10 
2.350xl0~ 10 
4.660xl(T 10 



t = 1.0 L x 
L 2 



3.936xl0~ 8 
5.048x10" 



3.80 
3.83 
3.85 



3.79 
3.83 



Loo 1.019xl0- 7 3.84 



3.90 
3.92 
3.93 



3.90 
3.92 
3.93 



t = 1.5 



Li 
L 2 

Loo 



2.400xl0~ 8 
3.078xl0~ 8 
6.212xl0~ 8 



3.79 
3.83 
3.84 



1.731xl0~ 9 
2.166xl0~ 9 
4.334xl0~ 9 



3.90 
3.92 
3.93 



1.158xl0- 10 
1.434xl0- 10 
2.842xl0~ 10 
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Tabic 6: Problem 3: L\-, L2-, L x - norm error in ui and spatial order of convergence with 
St = h 2 = k 2 



Time 



11 x 11 



Order 



21 x 21 



Order 



41 x 41 



£ = 0.5 Li 1.450xl0~ 6 3.85 
L 2 1.819xl(T 6 3.89 
Loo 3.480xl(T 6 3.87 



1.004xl0" 7 
1.226xl0- 7 
2.340xl0~ 7 



6.582xl0" 9 
7.946xl0~ 9 
1.501xl0~ 8 



3.93 
3.92 
3.96 



t = 1.0 



Li 

L 2 



;.900xl0~ 7 
1.116xl0- 6 



3.85 
3.89 



Loo 2.153xl0- 6 3.91 



6.172xl0~ 8 
7.530xl0- 8 
1.428xl0- 7 



3.93 
3.95 
3.96 



4.049xl0~ 9 
4.882xl0~ 9 
9.187xl0- 9 



t = 1.5 



Li 
L 2 

Loo 



5.423xl0~ 7 
6.795xl0" 7 
1.316xl0~ 6 



3.85 
3.89 
3.92 



3.763xl0- 8 
4.589xl0~ 8 
8.692xl0~ 8 



3.93 
3.95 
3.96 



2.469xl0~ 9 
2.976xl0~ 9 
5.593xl0~ 9 
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Table 7: Problem 4: L\- and L 2 - 
order of convergence 



norm difference in horizontal velocity and perceived 



Time Grid Spacing 
~25F — 



Lx-norm p L 2 -norm p 



h - k - 64 
St = 0.005 



1.0 h = k = & 
St = 0.0025 



U _ U _ 2lIL. 

ri — is, — 25g 
St = 0.0005 



=6.946xl0" 6 



\\U 3 -U 2 \\ 
=6.878xl0" 7 



3.19 



\\U3-U4 
=1.854xl0" 5 



\\U 3 -U 2 \\ 
= 1.092xl0" 6 



h — k — Mil 
11 - K - 64 

St = 0.005 



h — k — 
St = 0.0025 



4.00 



t = 5.0 



h = k 



2tt 
256 

at = 0,0005 

h = k = 
St = 0.005 



t = 10.0 h = k = ^ 
St = 0.0025 



ll^3-C/i|| 
= 1.766xl0" 3 



=1.354xl0" 4 



3.59 



11^3 "fill 

=3.203xl0" 3 



\\u 3 -u 2 \\ 

=2.992xl0" 4 



3.28 



11 — K - 256 

St = 0.0005 



\\Uz-u4 

=1.155xl0~ 2 



\\u 3 -u 2 \\ 

=2.498x10-3 



1.86 



ll^3-^i || 
= 1.627xl0" 2 



\\u 3 -u 2 \\ 

=3.962xl0" 3 



1.64 
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Table 8: Problem 4: L\- and L 2 
of convergence 



norm difference in vertical velocity and perceived order 



Time Grid Spacing 
~25F — 



Lx-norm p L 2 -norm p 



h - k - 64 
St = 0.005 



1.0 h = k = & 
St = 0.0025 



U _ U _ 2lIL. 

ri — is, — 25g 
St = 0.0005 



ll^s-K|| 
=2.041 xlO" 4 



\\V 3 -V 2 \\ 

:1.128xl0- 5 



4.10 



II V3 

=5.063xl0" 4 



||^3-^2|| 

=2.750xl0" 5 



h — k — Mil 
11 - K - 64 

St = 0.005 



h — k — 
St = 0.0025 



4.12 



t = 5.0 



h — k 



2tt 
256 

St = 0.0005 

^ = k = gf 
5t = 0.005 



t = 10.0 h = k = ^ 
St = 0.0025 



11^3 -1411 
=2.872 xlO" 3 



||^3-^|| 
=3.094xl0" 4 



3.05 



11^3 

=5.589xl0" 3 



||^3-^|| 
=8.369xl0" 4 



2.51 



11 — K - 256 

St = 0.0005 



11^3 

= 1.234xl0~ 2 



||^3-^2|| 

=2.305xl0" 3 



2.12 



11^3 

=1.733xl0" 2 



l|V 3 -V 2 || 
=3.494xl0" 3 



1.99 
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Table 9: Problem 5: Properties of Primary, Secondary and Tertiary Vortices for the lid- 
driven square cavity from Re = 1000 to Re — 7500. 



Vortex 


Property 


Ref. 


i?e =1000 


Re =3200 


Re =5000 


Re =7500 


Primary 






-0.1180 


-0.1214 


-0.1218 


-0.1220 






[29] 


-0.1179 


-0.1203 


-0.1190 


-0.1120 






[30] 


-0.1163 




-0.1142 


-0.1113 






[31] 


-0.117 


-0.122 


-0.122 


-0.122 




(x, y) 




0.5313, 0.5625 


0.5156, 0.5391 


0.5156, 0.5312 


0.5156, 0.5312 






[29] 


0.5313, 0.5625 


0.5165, 0.5469 


0.5117, 0.5352 


0.5117, 0.5322 






[30] 


0.5313, 0.5586 




0.5156, 0.5313 


0.5156, 0.5234 






[311 


0.5250, 0.5625 


0.5188, 0.5438 


0.5125, 0.5375 


0.5125, 0.5313 


Secondary 














TL 








8.10e-4 


1.60e-3 


2.36e-3 






[29] 




7.28e-4 


1.46e-3 


2.05e-3 






[30] 






1.75e-3 


3.14e-3 






[31] 




7.33e-4 


1.54e-3 


2.07e-3 




(x, y) 






0.0549, 0.8984 


0.0625, 0.9062 


0.0703, 0.9063 






[29] 




0.0547, 0.8984 


0.0625, 0.9102 


0.0664, 0.9141 






[30] 







0.0625, 0.9102 


0.0664, 0.9141 






[31] 




0.0563, 0.9000 


0.0688, 0.9125 


0.0688, 0.9125 


BL 


ib 




2.53e-4 


1.13e-3 


1.42e-3 


1.63e-3 






[29] 


2.31c-4 


9.78e-4 


1.36e-3 


1.47e-3 






[30] 


3.25e-4 




2.22e-3 


1.76e-3 






[31J 


2.02e-4 


1.03e-3 


1.32e-3 


1.60e-3 




(x, y) 




0.0782, 0.0781 


0.0804, 0.1203 


0.0703, 0.1406 


0.0625, 0.1563 






[29] 


0.0859, 0.0781 


0.0859, 0.1094 


0.0703, 0.1367 


0.0645, 0.1504 






[30] 


0.0859, 0.0820 





0.0664, 0.1484 


0.0703, 0.1289 






[311 


0.0875, 0.0750 


0.0813, 0.1188 


0.0750, 0.1313 


0.0688, 0.1500 


BR 






1.88e-3 


2.93e-3 


3.25e-3 


3.49e-3 






[29] 


1.75e-3 


3.14e-3 


3.08e-3 


3.28e-3 






[30] 


1.91e-3 




4.65e-3 


8.32e-3 






[3J 


1.70e-3 


2.86e-3 


2.96e-3 


3.05e-3 




0,2/) 




0.8594, 0.1094 


0.8203, 0.0859 


0.7969, 0.0703 


0.7813, 0.0625 






[29] 


0.8594, 0.1094 


0.8125, 0.0859 


0.8056, 0.0742 


0.7813, 0.0625 






[30] 


0.8711, 0.1094 




0.8301, 0.0703 


0.8828, 0.0820 






[31] 


0.8625, 0.1125 


0.8125, 0.0875 


0.8000, 0.0750 


0.7813, 0.0625 


Tertiary 














BR 


Ymin 






-3.90e-7 


-3.04e-6 


-6.35e-5 






[29] 




-2.52e-7 


-1.43e-6 


-3.28e-5 






[30] 






-2.47e-5 








[31] 




-2.37e-7 


-1.70e-6 


-1.89e-5 




(x,y) 






0.9880, 0.0114 


0.9744, 0.0211 


0.9453, 0.0469 






[29] 




0.9844, 0.0078 


0.9805, 0.0195 


0.9492, 0.0430 






[30] 






0.9668, 0.0293 








[21] 




0.9875, 0.0125 


0.9750, 0.0188 


0.9500, 0.0375 


Grid Size 




[29] 
[30] 
[31] 


65 x 65 
129 x 12§° 
257 x 257 

81 x 81 


129 x 129 
129 x 129 
257 x 257 
161 x 161 


129 x 129 
257 x 257 
257 x 257 
161 x 161 


129 x 129 
257 x 257 
257 x 257 
161 x 161 
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Figure 2: Computational stencil for the scheme: the used nodes are denoted by 
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Figure 4: Problem 2: Time evolution of errors corresponding to three different norms. 
Pe = 2,a = 100, c = d = 80, St = ti 2 = k 2 . 
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Figure 5: Problem 2: (a) Present, (b) PDE ^3/, (c) HOC 0/, (d) RHOC Pe = 200, 
a = 100, c = d = 10 4 , St = 2.5 x 10^ 5 , h= k = 0.02. Dash-dot contour lines in (a)-(d) 
correspond to the exact solution. 
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Figure 6: Problem 4: Gird convergence for velocity components at time t = 6.0 (denoted 
by dotted lines) and t = 10.0 (denoted by\$plid lines) (a) horizontal velocity along the 
vertical centreline and (b) vertical velocity along the horizontal centreline. 
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(b) 

Figure 8: Problem 5: Comparisons of computed profiles (denoted by solid lines) of (a) 
horizontal velocity along the vertical centreing: and (b) vertical velocity along the horizon- 
tal centreline with those of Ghia et al. pfij/ (denoted using symbols) for different Reynolds 
number. 
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